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Abstract 

Interval jitter and spike resampling methods are used to analyze 
the time scale on which temporal correlations occur. They allow 
the computation of jitter corrected cross correlograms and the per¬ 
formance of an associated statistically robust hypothesis test to de¬ 
cide whether observed correlations at a given time scale are signifi¬ 
cant. Currently used Monte Carlo methods approximate the proba¬ 
bility distribution of coincidences. They require generating Nmc sim¬ 
ulated spike trains of length T and calculating their correlation with 
another spike train up to lag r max . This is computationally costly 
0(Nmc x T x r max ) and it introduces errors in estimating the p value. 
Instead, we propose to compute the distribution in closed form, with 
a complexity of 0(C max log(C' max )r max ), where C max is the maximum 
possible number of coincidences. All results are then exact rather 
than approximate, and as a consequence, the p-values obtained are 
the theoretically best possible for the available data and test statis¬ 
tic. In addition, simulations with realistic parameters predict a speed 
increase over Monte Carlo methods of two orders of magnitude for 
hypothesis testing, and four orders of magnitude for computing the 
full jitter-corrected cross correlogram. 
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1 Introduction 


Though there are many means by which neurons communicate, using both 
chemical and electrical mechanisms, most attention has been paid to series 
of action potentials (spikes). It is known that in some cases the detailed 
time structure of these spike trains is used for information transmission 
while in others, only the overall number of spikes in some interval seems 
to be important but not their position in the interval. Examples of the first 
type are various kinds of “temporal coding” schemes proposed for different 


Singer and Gray 

(1995) 

Riehle et al. 

(19971); 

_ . 7 * _ ~ ^ 

Niebur et al. 

( 

19931): 

Softkv 

(1995 

1: Steinmetz et al. 

(2000), while the latter is the well-known rate-code 


(119981) . To distinguish between these two possibilities, it is necessary to find 


whether reproducible correlations at the relevant time scale are present in 
neuronal data. One common way to approach this problem is to use auto- 
or cross-correlation functions as test statistics. Then one can (a) search for 
non-trivial structure in the function, like deviations from uniformity, and 
(b) detect whether there are differences between these functions in different 
experimental (e.g. behavioral) conditions. 

The situation is complicated by the influence of rate variations on the raw 
correlations. To increase the signal-to-noise ratio, correlations are typically 
computed as averages over many trials. Changes in the behavioral state 
of an animal, e.g. due to onset of sensory stimuli or motor responses that 
occur always at the same time during a trial, typically result in changes in 
neural firing rates which are common to many neurons. While these are 
genuine correlations, they are unrelated to the neuronal coding question. 
Different techniques have been develo ped to remove them, e.g. subtraction 


of a “shuffle predictor” Perkel et al 


( 19671 ). the average of cross correlations 
between spike trains from permuted trialqj. While this correction removes 
correlations that are time-locked to trial onset, it was later pointed out that 
peaks in the correlation function that may be taken as indicative of correla ted 
firing (e.g. at zero lag) can also be caused by slow rate covariations [Brod y 
( 1998 . 1999() . After finding a significant peak in the cross correlation function, 
this ambiguity can be addressed by analyzing the time scale at which the 


1 A “shift predictor” is very similar but the correlation function is computed from trials 
that immediately follow each other, rather than from randomly selected trials. 
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measured correlation arises. 


It was pointed out more recently I Amarasingliam et al.l (12012) that the 


null hypothesis of spike trains being independent in earlier work Perkel et al 


( 19671 ) is useless if their mean rates co-vary. Rate co-variation of means the 
spike trains are not independent which leads to its immediate rejection of the 
null without providing any further insight. Amarasingliam and his colleagues 
instead proposed a more detailed null hypothesis, namely that within an in¬ 
terval of width A, the exact location of spikes does not matter. Then, under 
the null hypothesis, simulated spike trains can be generated by modifying 
the spike times of an original measured spike train within a range of A. The 
cross correlations obtained from these modified (“jittered”) spike trains are 
then compared to those obtained from the original. If significant differences 
are found, the null hypothesis is rejected and it is likely that non-random 
correlations at time scales < A are present in the data. Additionally, this 
method gives rise to the computation of jitter-corrected cross correlograms, 
which have be en used to compare change s in synchrony acro ss experimen¬ 
tal conditions ( Hirabavashi et al. . 2013a bl: Smith et al.l . 2013 ). Because the 
method relies on repeated simulation of spike trains, it will be referred to as 
the Monte Carlo jitter method for the purposes of this paper. 

While the Monte Carlo jitter method (described fully in Section 12.ip is 
useful and easily generalized to complex statistical tests and hypotheses, 
its practical utility is limited by the inherent trade-off between accuracy and 
computation time in all Monte Carlo methods. As we will show in Section 13721 
the computation time may be prohibitively long, and even at this cost, it will 
only be a numerical approximation of the true solution. In the case where 
the test statistic is the cross-corrcllation value at a single lag, the p value 
can be computed exactly, as was shown by Harris on 020131 ). In the present 
study, we therefore explore the benefits of computing in closed form the 
distribution which is only approximated by the Monte Carlo simulations. 
Accordingly we refer to this method, described in Section 12.21 as the closed 
form jitter method. In addition to computing the p value for rejecting the 
null hypothesis exactly we show that the computation of the jitter-corrected 
cross correlogram follows readily from that derivation. The computational 
performance of the closed form jitter and Monte Carlo jitter methods are 
compared theoretically (as computational complexity) in Section 13.11 and 
practically (as computational time) in Section 13.21 
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T 

X, Y 

C(t) 

A 


J 

N(X,j) 

YMC 

Xmc 

qMC 

R t 

Pt 

JCCG(r) 

P t (C mc ) 

Pt(C(t)) 

7~max 

X 

a 


max 

max 


Length of binned spike trains 
Two spike trains 
Correlation Lag 
Correlation of A" and Y 
Jitter interval width 

Index over Monte Carlo simulated signals 
Index over jitter intervals 
Number of spikes in X in interval j 
ith Monte Carlo simulated signal 
Number of Monte Carlo simulated signals 
Correlation of A, MC and Y 
Number of cases where C^ c (r) > C(r) 
p Value for correlation at lag r 
Jitter-corrected Cross Correlogram 

Monte Carlo estimate of the distribution of correlations at lag r 

Number of coincidences in the jth interval 

True distribution of correlations at lag r 

Maximum r value processed 

Maximum value of N(X,j) or N(Y,j) 

Maximum possible number of coincidences 


Table 1: Glossary. Variables are listed in the order in which they are 
introduced. 
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2 Materials and Methods 


2.1 The Monte Carlo Jitter Method 


Utilizing the Monte Carlo jitter method flAmarasingham et all 2012), it is 
possible to determine whether correlations arise from fine temporal structure 
or larger scale variations, sometimes referred to as rate covariations. This 
determination is made by comparing a test statistic (in this case cross cor¬ 
relations) of an original pair of spike trains against those computed from a 
set of jittered spike trains as described below. The jitter method, like cross 
correlation, operates on binned spike trains which we take as binary signals 
with values 0 and 1 and integer arguments 0 to T — 1, where T is the number 
of bins in the spike train. The binary assumption implies that the bin size is 
small enough (typically 1ms or so) such that two spikes cannot be recorded 
in a single time bin. A sufficiently small bin size can always be chosen since 
there are limits on the minimal inter-spike interval time due to the refractory 
period of the neurons in question. 

Let X(t ) and Y(t ) be two such binned spike trains. The processing then 
consists of the following steps: 


1. Compute the cross correlation C(r) between the original X and Y, 

c(T) = £x(«-T)y(t) 

t 

where the sum runs from 0 to T — 1 and X is assumed to be 0 if its 
argument is outside that range. 

2. Subdivide one of the signals, say X, into intervals of width A. 

3. Count the number of spikes in each interval of X. In interval j this is, 

Cj+1)A—1 

N(x,j)= Y. (1) 

k=jA 


4. For X generate Nmc Monte Carlo simulated signals {A" 4 MC }, in which 
the spike counts for each interval are the same as in the corresponding 
interval in A", such that N( A 4 MC ,j) = N(X,j ) for all i,j. However, 
now spike times within the interval are all equally likely. Spike times 
should be sampled without replacement to ensure that the spike count 
stays constant without putting multiple spikes in a single bin. 
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5. Compute the cross correlation C' l MC (r) for lag time r between each 
Xj MC and the second spike train Y to get an estimate of the distribution 
P t (C mc ) of cross correlation values for each time lag r. 

C, mc (t) = ^A', mc («-t)K(*) 

t 

6. Let R r be the number of simulations where C J MG (r) > C'(r). Then the 
p value for a given lag r is computed as 

R t + 1 
^ Nuc + 1 

7. If desired, a jitter-corrected cross correlogram (JCCG), defined using 
the expectation operator E[-], can be computed as 

JCCG(r) = C(r) - E[C M0 (r)] « C(r) « A- £ C“°(r) (2) 

™MC 


where the approximation approaches equality with Nmc 00 • 

Jittering the spikes within an interval of size A destroys all correlations 
at time scales within this interval. The cross correlations computed from the 
jittered spike trains therefore are not correlated on time scales A or smaller, 
and P t (C mc ) is the distribution of correlations at time lag r obtained under 
the null hypothesis that correlations at time scales < A are indistinguishable 
from random correlations. If the measured cross correlation C(r) is signif¬ 
icantly outside this distribution, we have to reject the null hypothesis and 
we conclude that nonrandom correlations at lag r with time scales < A are 
found in the observed spike trains. If, on the other hand, the observed cor¬ 
relation is consistent with what is seen in the distribution of jittered spike 
trains, then we cannot reject the null hypothesis. This means we cannot 
exclude that the observed synchrony is caused by correlations on time scales 
outside the jittered range, in other words that the observed correlation at lag 
time r is caused by rate variations on time scales greater than A. 

In practice, X(t) and Y(t) do not have to be gathered from a continuous 
block of time. In the case of multiple trials of the same experimental con¬ 
dition, it may be useful to concatenate the recorded spike trains (possibly 
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after removing sections of them, like those recorded during stimulus onsets). 
In doing so, a period of no spiking of width r max (the largest correlation lag 
of interest) should be added between the trials so that correlations between 
trials don’t affect the outcome of the jitter procedure. 

As mentioned in Section [Q the practical utility of the Monte Carlo method 
is limited by the trade-off between accuracy and computation time inherent 
in all Monte Carlo algorithms. Furthermore, in practice a single set of Monte 
Carlo simulations is often generated for many hypothesis tests (he. tests at 
multiple lags), introducing potential dependencies between the different tests 
when they should be treated independently In order to avoid both of these 
issues, the probability distribution P r (C MC ) can be computed exactly and 
independently for each time lag as described in the following. 

2.2 Closed Form Computation 

2.2.1 Probability Distribution For One Interval 

First, let us consider a single interval consisting of A time bins. For example, 
if spike times have been binned to 2 ms, for an interval of width 20 ms we 
have A = 10. Since time has been discretized, it is still possible to discuss 
this unitless value as a length of time, a time scale, or a interval width for 
a given bin size. As before, we assume that the sequence is binary, so each 
bin has either zero or one spike. This is true even when the spike times 
are jittered because spike times are sampled without replacement. For this 
single interval, the probability of a given number of coincidences occurring is 
determined by three values: A, N(X,j) the number of spikes in interval j of 
spike train A", and N{Y,j) the number of spikes in interval j of spike train 
Y. 

As a first step, we count the number of perfect coincidences, in which one 
spike occurs in both X and Y within the same time bin, meaning r = 0. 
Using the standard notation of (“) for the combinatorial operation (a choose 
6), we find that there are ( N ^^) ways to distribute N(X,j ) spikes in A 
available bins. The number of empty (spike-less) bins in spike train Y is 
[A — N(Y,j)]. The number of ways to distribute N(X,j ) spikes such that 
each of them falls into one of these empty bins is (^pryj 7 ^)- These are all 

2 The procedure of generating one set of spike trains for multiple lags is appropriate 
inappropriate only if each lag is being tested independently. If the test statistic is the sum 
of C(r) over a range of lag values, a single set of simulated spike trains is appropriated 
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possible cases in which a coincidence is avoided. The probability that zero 
coincidences occur in the j'-th interval is therefore 


P(C%* = 0\A,N(X,j),N(Y,j)) 


1 N(X,j) ) 


(3) 


where Cj nt is the number of coincidences in this interval. 

We can generalize equation [3] to a non-zero number c of coincidences 
by breaking the numerator up into the number of ways that c spikes can 
coincide with the spikes in Y, and N(X,j) — c spikes coincide with the gaps 
(or non-spikes) in Y. We can thus compute a probability distribution for 
each interval j, 


/A -N(Y,j)\ /N(Y,j)\ 

P(Cf = c|A, N(X,j), N(Y,j)) = (A ', a\‘ W 

(n ( xJ 

where we follow the customary convention of setting the value of a “choose” 
operation to zero if either of its arguments is negative, or if its upper argu¬ 
ment is less than the lower. If this happens in the numerator of eq. [4l the 
probability on the left hand side becomes zero. Of course, the denominator 
is always positive since N(X,j) < A. This is a hypergeometric distribution. 

Equation [4] is easily generalized to nonzero values of r by applying the 
analysis leading to it to a shifted version of Y. For the computation of 
N(Y,j ), this implies adding r to the summation limits in eq. CD As with 
other correlation algorithms, the boundaries of finite spike trains (beginning 
and end) result in fewer intervals to analyze as r gets further away from 
zero. Thus generalizing eq. [I] to non-zero r, we denote the resulting number 
of coincidences as C‘ nt (r) and the associated probability distributions as P) nt . 


2.2.2 Jitter-Corrected Cross Correlation 


Once we have the analytical probability distribution for the correlations, we 
can obtain all relevant quantities to characterize the pairwise correlations be¬ 
tween two spike trains. It is straightfo rward to compute the commonly used 


jitter -corrected cross correlogram (e.g.. lHirabavashi et al.l.l2013alJbl: ISmith et al 


20131 ) which shows the correlation function after all correlations on time scales 
longer than A have been removed. It is defined in analogy to equation [2] 
where the expectation value of the stochastic solution, P[C MC (r)], was used. 
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We can replace this approximation by the exact solution A[C' mt (r)]. Fur¬ 
thermore, by the null hypothesis each interval is conditionally independent 
based on the spike counts. Therefore, the JCCG can be computed without 
approximation by 


JCCG(r) 


C(t) — E 


£c«(t) 


L J 


j 


(5) 


which, as should be remembered, is computed for a specific jitter interval 
width A. The expectation on the right can either be calculated for each 
window as N(x,j ) x N(y,j)/A. 

The jitter-corrected cross correlogram is used, for instance, when the 
scientific question of interest is whether there are significant changes in syn¬ 
chrony between conditions, rather than a test of the presence or absence of 
synchrony. It is then used as part of a bootstrap statistical test in which 
the observed pairwise correlation is compared with the distribution obtained 
from eq. 0 


2.2.3 Probability Distribution For Spike Train 

One can also obtain the probability distribution for the entire signal P t (C(t )) 
as the convolution of the individual probability distributions for all intervals, 
?j nt . This is identical to computing P r (C MC ) from Section I27T1 with an infinite 
number of Monte Carlo simulations for each value of r. One can then evaluate 
how likely it is that the observed cross correlation C{r) is explained by this 
probability distribution. The likelihood p that this is the case is obtained as 
the integral of the probability density function exceeding G(r), as in 

OO 

Pr= Pr ( C ) ( 6 ) 

c=C(r) 


3 Results 

3.1 Computational Complexity 

In many situations, the statistical distributions underlying the phenomena 
under study are complicated or unknown and performing Monte Carlo simu¬ 
lations are the only way to make progress, even though it may be costly and 
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it introduces additional randomness in the processing. In the case considered 
here (binary spike trains, null hypothesis of uniform spike time distribution 
in fixed interval, cross correlation test statistic), the distribution P T (C ) can 
be computed directly, using the closed form jitter method described above, 
without the need for repeated simulations. This section will compare the 
computational complexity of using the Monte Carlo jitter method against 
the direct computation using the closed form jitter method. 


3.1.1 Monte Carlo Method 


In the Monte Carlo algorithm, the data generation step requires a permuta¬ 
tion of A data points for each interval and simulation. Since a single per¬ 
mutation operation has a computational complexity 0(A), and A times the 
number of intervals is the length of the signal T, generating the set of Monte 
Carlo simulations {A 4 MC } is 0(Nmc x T). The complexity of cross correlation 
or convolution of two signals with lengths T is 0(T x logT), assuming an 
FFT-based method f Coolev and Tukev . 19651 ) is used. So computing the full 
Monte Carlo probability distribution for all values of r is 0(N MC xTx log T ). 
In many cases, not all values of r are needed. If the correlation is computed 
only for the subset of delay values from 0 to r max , the complexity for the 
Monte Carlo jitter method is 


0(Nmc x T x T max ) 

Computing the jitter-corrected cross correlogram by this method only 
requires one additional sum, with complexity 0(Nmc x r max ) so the total 
complexity remains unchanged. 


3.1.2 Closed Form Probability Distribution 

To compute the exact probability distribution with the closed form jitter 
method, note that the values of the distribution can be precomputed based 
on the maximum values of N(X,j ) and N(Y,j) over all j; call this maxi¬ 
mum JVjpax. Also, n choose k operations can be as fast as 0(min(A;,n — k )) 
( Manolopoulosl . .2002). Therefore, a three dimensional table of all possible 


values of P(Cj nt \N(X, j), N(Y, j)) can be precomputed and then looked up 
for each interval. Generating this table requires up to A ni ay different values 
of C, iV max values of N(X, j), and iV max values of N(Y,j). Computing each 
value requires three n choose k operations, which are on the order of 0(N max ), 


11 












Closed Form Jitter Analysis 


so the total computation of the probability table is 0(N m ax 4 ). While the ex¬ 
ponent is high, the expression does not have any dependence on the length 
of the signal and, furthermore, iV max < A is a small number in essentially 
all cases of interest. In practice, for analyzing neurophysiological data it is 
rare that a time resolution finer than 1 ms is needed, or controlling for cross 
correlations at time scales larger than approximately 100 ms (i.e. A ss 100). 
The full lookup table is therefore maximally a 100 x 100 x 100 matrix, which 
requires negligible resources to compute and store. 

To compute the combined probability distribution P T (C ) over all inter¬ 
vals, all interval probability distributions P r (Cj nt ) must be convolved, and 
the computational complexity of the problem is dominated by these convo¬ 
lution operations. As will be shown, we can improve performance by taking 
advantage of the structure of the problem at hand, since many of the convo¬ 
lution operations are identical. As a result, the convolutions can be grouped 
together based on N(X,j ) and N(Y,j) and quickly combined so that 0{T ) 
convolutions will turn into 0(-/V max 2 3 4 5 ) convolutions. This can be done by the 
following procedure: 

1. Take the Fast Fourier Transform (FFT) of P T (CJ lt \N(X, j), N(Y, j)) 
for each encountered value of N(X,j ) and N(Y,j). 

2. Raise each complex frequency spectrum value to a power equal to the 
number of times that the (N(X, j), N(Y, j)) pair appears. 

3. Multiply these frequency spectra. 

4. Take the inverse FFT of the result to get the final probability distribu¬ 
tion P T (C) and compute p T as in equation [6l 

5. Repeat steps [2] through H] for each value of r to be tested. 

The FFT operations in step [T] must be zero-padded up to the maximum 
number of coincident spikes C m ax to account for the highest possible number 
of synchronous spikes in the combined probability distribution. Therefore the 
FFT operation in step [1] is 0(C max x log(O max )). In step [2] exponentiation is 
0(1). However there are 0(T x A max 2 ) exponents to be taken, repeated r max 
times in step [3 Step |3] requires 0(T x A max 2 ) multiplications, again repeated 
An ax times. In step [2 the length of the spectral signal (to be inverted by FFT) 
is C max , so the operation is 0(O max x log(C max )) repeated r max times. When 
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combining these steps, note that C ma y is proportional to T. However C max 
will be used when relevant because it captures the frequency dependence of 
the computation time. Therefore the total computational complexity is 

C(C max x log(C max ) x T max ) 

Note that because the zero-frequency component of a probability distribu¬ 
tion is always exactly unity, the inverse FFT computation will have accuracy 
limited by the precision of the numerical system. In practice this implies 
that p values less than 10~ 13 will not be estimated accurately. 

3.1.3 Jitter-Corrected Cross Correlation 

Both the complexity analysis and the actual computation of the jitter-corrected 
cross correlogram is much simpler than that of the probability distribution. 
We generate a lookup table of possible E [C 111 *] values and, from equation [5J 
the jitter-corrected correlogram can be computed at a speed of 

0(T x r max ) 


3.2 Computational Execution Time 

For practical applications, consumption of resources is an important limita¬ 
tion for any computational method. For the size of problems encountered in 
typical neurophysiological experiments, the only limiting resource is execu¬ 
tion time. To compare the performance of the Monte Carlo jitter and the 
closed form method, the two algorithms were run side by side in the MAT- 
LAB environment (Mathworks, Natick MA). Synthetic spike trains were gen¬ 
erated that varied in both frequency of spiking (5 to 500 Hz) and length (1 
to 91 seconds). For each (time, frequency) condition, 50 spike trains were 
generated, binned to 1 ms, and the average processing time was computed. 
Processing was performed with r max = 100/ns and A = 20. All computations 
were performed on an Intel i7 920 processor with 12 GB of RAM running 
Linux Ubuntu 12.04. 

For the Monte-Carlo method, Amc was set to 1000. This selection of 
Amc is unrealistically low for two reasons. First, it can at best result in a 
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Bonferroni corrected p value of 0.201 due to the 201 p values being tested 
in the range of — r max to r max . As the execution for A MC = 1000 already 
takes 5.7 days to run. increasing Amc is impractical. Second, only a single 
set of Monte-Carlo trials were generated for all lag values computed, induc¬ 
ing potential correlations between the p values. These correlations should 
decrease as more trials are generated. Therefore results are extrapolated to 
Amc = 20, 000 (resulting in a minimum p ss 0.01) under the assumption 
that the processing for 20 times as many simulations would take 20 times as 
long. Though the bonferroni correction used here is conservative, it is less 
conservative (by an order of magnitude) than simulating a whole new set of 
spike trains for each p value as would be required to entirely eliminate any 
correlations between the p values. 

For the closed form jitter method, all lookup tables were computed de 
novo for each spike train. This is a conservative approach (favoring the 
Monte Carlo technique) since performance of the closed form jitter method 
could be improved by computing the tables only once and using them for all 
spike trains. This is certainly advised in a “production environment.” 

The results of this simulation, shown in Figure [Q illustrate a number of 
features about the speed of the two algorithms. Plotted is the performance 
gain, defined as the ratio of the computation time between the Monte Carlo 
jitter method and the closed form jitter method. The first observation is 
that the closed form method is substantially faster than the Monte Carlo 
method in all cases considered. Second, while the performance gain depends 
only weakly on spike train length, it does decrease with increasing bring 
rate. This is because the computation time of the closed form jitter method 
increases with bring rate. In practice, however, it is rare to observe bring 
at sustained frequencies exceeding 100 Hz in physiological recordings. In the 
physiological range, the closed form jitter algorithm is faster by a factor of 
approximately 180 to 7200. 


F K 

Har rison (120131 ) uses importance sampling to accelerate the Monte Carlo 
hypothesis testing process which requires drawing fewer samples. In that 
work the number of samples needed, even for a low Bonferroni corrected p 
value, is reduced to 100. However, each sample is reported to take 18 times 
as long to generate and process as before, ebectively resulting in a simulation 
about 11 times faster than the Monte Carlo simulation with Amc — 20, 000. 
Therefore, under physiological conditions the closed form computation has an 
expected speed-up of 16 to 650 times compared to the importance sampling 
method. It should be noted that in cases where even lower p values are needed 
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because of multiple hypothesis constraints, importance sampling will provide 
larger gains in estimating very small p values. In these cases, increasing the 
p value requirements has no effect on the computation speed of the closed 
form method, so the closed form method can be expected to be faster in all 
cases. 

Another improvement mentioned in Section 12.2.21 is the ability of the 
closed form jitter method to compute the jitter corrected correlogram very 
rapidly, without computing the null hypothesis distribution of correlation 
values. To show the magnitude of the improvement, the simulation was 
repeated with only the mean of the null hypothesis distribution computed 
under the closed form jitter method since this is all that is needed for the 
corrected correlation function, eq. [5] We also restricted firing frequencies to 
the range 5-200 Hz. Figure [2] shows the ratio of the time it takes to compute 
equation [2] vs. equation [5] In these cases, the closed form jitter calculation 
is substantially faster (480x-13,000x), with increasing benefits for increasing 
spike train lengths. As discussed previously, the spike train length is typically 
not that of individual trials but of the concatenation of many trials. 


4 Discussion 


The importance, or absence of it, of precise timing of neural spikes has been 
discussed for the last half-century. Several techniques have been developed 
to characterize neuronal responses at fine time scales and it is clear that 


Gawne and Richmond! 

1993; 

Rov et al.. 

2000) 


conclusions (e.g. 

tant difficulty is that firing rates can co-vary in the neurons under study. It 
is well-known that such co-variations are observable in quantities like pair¬ 
wise cross correlation functions but they are typically considered as irrelevant 
from the point of view of neuronal coding or of determining the connectivity 
in the underlying circuitry. For instance, the onset of a stimulus will typi¬ 
cally generate a temporary increase in firing rates in sensory cortex but the 
resulting increase in cross correlation is usually n ot cons id ered o f importance 
for neural coding (for an exception see Chase and Young . 20071 who showed 
that spike timing relative to onset-related population activity is informative). 
One common way to subt ract su ch stimulus-locked effects is by subtracting 
a “shuffle predictor” ( Perkel et ah . 1967t h obtained by computing cross cor¬ 
relations between spike trains from permuted trials. It has been pointed out 
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repeatedly (see references in the Introduction) that this does not eliminate 
spu rious correlations, inclu ding close to r = 0 (synchron y). 

Brodv f 19981. 19991) and Amarasingham et al. ( 2012 1 proved that adopt¬ 


ing the null hypothesis of independent neurons can not solve the problem. 
Observation of such correlations is, indeed, evidence against the null hypoth¬ 
esis of independence between the two observed spike trains. Rejection of this 
null hypothesis can, however, occur either because spikes in the two spike 
trains are correlated “one-by-one” (synchrony), or because of slow bring rate 
covariations common to both spike trains. The fact that this null hypothesis 
can be rejected does not tell us why it is rejected. If the question is whether 
synchrony exists at less than a given time scale (only), this is the wrong null 
hypothesis. Instead, the time scale n eeds to b e sp ecibed explicitly. The null 
hypothesis chosen by Am arasingham et ah (2012) is that changes of spike 
times within a time interval of size A have no ebect on the computed statis¬ 
tic, in this case the correlation function. It_is_this null hypothesis that is 
tested by computer simulation in the Amarasingham et ah ( 2012f) study and 
analytically in this report. 

A key element of the methods discussed here is that the jitter intervals 
are dehned without reference to the original spike trains. This ensures that 
if the null hypothesis is true, there is no way to distinguish the original 
spike trains from the Monte Carlo simulated spike trains. This characteristic 
(called exchangeability) ensures that the obtained p values are from a well 
formulated hypothesis test. If, on the other hand, the resampling method 
was changed so that each spike was jittered about it’s original spike time, 
then even under the null hypothesis the original spike train would stand out 
from the rest because all of its spikes would be at the center of the jitter 
intervals. Therefore the resulting test w ould nqt^be a proper statistical test 
and should be avoided (1 A mara sing ham et ah, 2012). 

We have discussed two ways one can choose to characterize the correla¬ 
tions between two spike trains. One is a strict hypothesis testing approach. 
A null hypothesis is formulated, namely that the observed correlations are 
indistinguishable from correlations between spike trains whose spikes have 
been distributed randomly within intervals of length A, without changing 
the number of spikes in each interval. By comparing the observed correla¬ 
tion with those in the distribution generated under the null hypothesis, it is 
then decided for a given a whether the null hypothesis can be rejected. 

The alternative is to compute the time-resolved correlation function and 
“correct” for the correlations as observed under the null hypothesis, by sub- 
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tracting the expectation value of the latter. This is the more commonly cho¬ 
sen approach, perhaps because the time-resolved correlation function is both 
intuitive and familiar. The distribution of JCCG values can be compared 
between experimental conditions (indicating a change in ’excess synchrony’) 
using a bootstrap test to test for significance. Also, its shape (e.g. the lo¬ 
cation of peaks) may provide insight that goes beyond the yes-no answer 
whether the null hypothesis can be rejected or not. 

In the Amarashigharri et ahl (2012) study, the Monte Carlo procedure is 
further developed to account for more potential causes of fine timing effects 
besides synchrony such as ramping spike rates within an interval or inter¬ 
spike interval distribution effects. These methods are straightforward and 
statistically well-defined. Like any Monte Carlo method, however, they only 
generate an approximation to the underlying distribution whose quality de¬ 
pends on the number of surrogate spike trains. In practice what is more prob¬ 
lematic is that the method can be computationally very costly. For instance, 
as discussed in section 13.21 our example problem using the simplest of the 
null hypotheses discussed (50 spike trains of a few seconds long each, mean 
rates between 1 and 100 Hz, maximal time lag of 100 ms, a = 0.01 with Bon- 
ferroni correction applied) would have required a simulation several months 
long on a reasonably fast machine. We therefore only simulated Nmc — 1000 
trials and extrapolated to the execution time needed for Nmc = 20, 000 but 
even that abbreviated Monte Carlo run took nearly six days. Some progress 
can be made by using much faster machines or many machines (the problem 
parallelizes easily) but execution time is clearly a problem. 

In contrast, the closed form jitter methods this report focuses on are ex¬ 
act, rather than approximate. More important for practical applications may 
be that they are extremely efficient, with a speed-up of at least two orders of 
magnitude for the hypothesis testing approach, and four orders of magnitude 
for_the full corr elation functions. Even over importance sampling methods 
( Harrison . 20131) . they have been shown to provide a substantial increase in 
speed. For the hypothesis testing examples used in our study (whose scope is 
quite comparable to that of typical neurophysiological experiments, assum¬ 
ing a proper Bonferroni correction is applied), computation time is reduced 
from more than 100 days under the original Monte Carlo method to about 
one night. Computational time required for the full correlation function is 
reduced from over 100 days to a few minutes. An increase in performance on 
this scale is more than merely a quantitative improvement. For instance, it is 
essentially impossible to explore variations in the analyses (like the influence 
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of the jitter time scale A) if each computational run takes a few months, but 
it is easy to do if it takes minutes. 

So far we were only concerned with correlations between two spike trains. 
Modern recording techniques are already increasing the number of simulta¬ 
neously recorded spike trains to tens or hundreds. Unfortunately, the closed- 
form jitter method is limited in the ability to analyze large ensembles. This 
is because the correlation functions of some pairs in an ensemble will restrict 
the possible correlation values of other pairs. For example, if there are three 
neurons X, Y, and Z, and the pairs XY and YZ have perfect correlation, 
then the pair XZ must also have perfect correlation. A Monte-Carlo jitter 
analysis that jitters an entire ensemble of neurons and then performs a hy¬ 
pothesis test on the ensemble can be performed relatively simply, but no such 
closed-form method exists yet. In order to avoid the constraints of the type 
described above, only N — 1 pairs of neurons can be analyzed with closed 
form methods when N neurons are recorded. 

Additionally, the nature of the exact solutions provides an opportunity 
for further exact analysis. Having a closed form solution allows questions 
about the effects of spike sorting errors, the value of A, or the structure of 
JCCG(t) to be addressed rigorously and more precisely than is possible with 
any numerical method. 

In conclusion, we study a statistical framework for quantifying correla¬ 
tions between spike trains at given time scales. It can be applied both for 
hypothesis testing and for correcting observed correlation functions for cor¬ 
relations at these time scales. Results are exact, and both computational 
complexity and computational time for realistic examples are several orders 
of magnitude lower than related approaches based on Monte Carlo simula¬ 
tions. 

Matlab code will be made available by the authors upon request. 
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Figure 1: Performance gain in implementing the closed form jitter method 
for p value computations. Gain is defined as the ratio in computation time 
between the Monte Carlo Jitter method and the closed form jitter method. 
Processing parameters used are r max = 100ms, A = 20, and lV M c = 20, 000 
(see text for details). 
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JCCG Computation 



Figure 2: Performance gain in implementing the closed form jitter method 
for Jitter Corrected Correlogram computations. Gain is defined as the ratio 
in computation time between the Monte Carlo Jitter method and the closed 
form jitter method. Processing parameters used are r max = 100ms, A = 20, 
and Nmc = 20,000 (see text for details). The lowered performance gain at 
signal length of 1 second is due to the overhead of computing the probability 
table de novo for each spike train. 
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